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A general method to treat non-Gaussian vacuum wave functionals in the Hamiltonian formulation 
of a quantum field theory is presented. By means of Dyson-Schwinger techniques, the static Green 
functions are expressed in terms of the kernels arising in the Taylor expansion of the exponent of the 
vacuum wave functional. These kernels are then determined by minimizing the vacuum expectation 
value of the Hamiltonian. The method is applied to Yang-Mills theory in Coulomb gauge, using 
a vacuum wave functional whose exponent contains up to quartic terms in the gauge field. An 
estimate of the cubic and quartic interaction kernels is given using as input the gluon and ghost 

■ propagators found with a Gaussian wave functional. 
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jy*} . I. INTRODUCTION 

(T) ' . 

£Sj ■ According to our present understanding of nature, Quantum Chromodynamics (QCD) is the theory of the strong 
interaction. This theory has been tested in the high-momentum or ultraviolet (UV) regime, where perturbation 
theory is applicable due to asymptotic freedom. Our knowledge on the low energy, strongly interacting regime of 
QCD stems mainly from lattice calculations, which have at least qualitatively reproduced many physical observables, 
in particular, the linearly rising confining potential for heavy quarks. Furthermore, lattice calculations have revealed 
the relevance of topological field configurations such as magnetic monopoles and center vortices for infrared phenomena 
like confinement and spontaneous chiral symmetry breaking. These calculations support the dual Meissner effect and 
> the vortex condensation picture of confinement [1]. Despite these substantial physical insights provided by the lattice 
I \ calculations, a thorough understanding of these infrared phenomena will not come from lattice calculations alone but 
j> . will require also studies of the continuum theory. 

Q\ ' In recent years there have been substantial efforts devoted to a non-perturbative treatment of continuum Yang-Mills 
Q\ , theory. Among these are a variational solution of the Yang-Mills Schrodinger equation in Coulomb gauge [2-5]. In 
this approach, using Gaussian type wave functionals, minimization of the energy density results in the so-called gap 
equation for the gluon energy (or static gluon propagator) . This equation has been solved analytically in the infrared 
^\ • [6] and in the ultraviolet [7] and numerically in the full momentum regime [4, 8]. One finds a gluon energy, which in 
the UV behaves like the photon energy but diverges in the infrared (IR), signalling confinement. The obtained gluon 
energy also compares favourably with the lattice data [9]. In particular, the infrared regime is correctly reproduced, as 
far as we can tell from available lattice data. There are, however, deviations in the mid-momentum regime (and minor 
. ones in the UV) which can be attributed to the missing gluon loop, which escapes the Gaussian wave functionals. 
These deviations are presumably irrelevant for the confinement properties, which are dominated by the ghost loop 
(which is fully included under the Gaussian ansatz), but are believed to be important for a correct description of 

■ spontaneous breaking of chiral symmetry [10]. 

The numerical wave functional obtained from the variational solution of Ref. [8] seems to embody the correct 
infrared physics as is revealed in the various applications considered to date: one finds a linearly rising static quark 
potential [8] , an infrared enhanced running coupling constant with no Landau pole [6] , a topological susceptibility in 
accord with lattice data [11], a perimeter law for the 't Hooft loop [12], and, within an approximate Dyson-Schwinger 
equation, an area law for the spatial Wilson loop [13]. Furthermore, in Ref. [14] it was shown that the inverse ghost 
form factor of Coulomb gauge Yang-Mills theory represents the dielectric function of the Yang-Mills vacuum and 
the so-called horizon condition [15] (of an infrared diverging ghost form factor) implies that the Yang-Mills vacuum 
is a perfect color dielectricum, i.e., a dual superconductor, which establishes the connection between the Gribov- 
Zwanziger confinement scenario [15, 16] and the monopole condensation picture [17, 18]. Finally, in Ref. [19] the 
functional renormalization group flow equation of the Hamiltonian approach to Coulomb gauge Yang-Mills theory 
was studied, yielding results for the gluon and ghost propagator similar to that of the variational approach [4]. 

In the present paper, we generalize the variational approach to the Hamiltonian formulation of Yang-Mills theory to 
non-Gaussian wave functionals. We will present a general method to treat non-Gaussian wave functionals in quantum 
field theory. The method is based on the observation that expectation values in the Hamiltonian formulation of 
d = 3 + 1 dimensional quantum field theory can be formally obtained from a generating functional of d — 3 dimensional 
Euclidean quantum field theory with an action defined by the logarithm of the vacuum wave functional. Expanding this 
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action functional in powers of the underlying field results in "bare" n-point kernels 7 n as expansion "coefficients" . We 
then exploit Dyson-Schwingcr equation (DSE) techniques [20-22] to express the expectation value of the Hamiltonian 
(if) in terms of these kernels j n , which are then determined by the variational principle, i.e., by minimizing (if). This 
approach is then applied to the Hamiltonian formulation of Yang-Mills theory in Coulomb gauge to include three- 
and four-gluon interaction kernels in the exponent of the Yang-Mills vacuum wave functional. 

By using such a non-Gaussian wave functional, the gluon loop is retained in the expectation value of the Hamiltonian. 
Although the gluon loop is irrelevant for the IR properties, it certainly influences the mid-momentum and UV regime 
of the gluon propagator and thus of the running coupling, and also contributes to the anomalous dimensions. As a first 
estimate of the effects of the non-Gaussian terms in the wave functional, we will calculate the gluon-loop contribution 
to the gluon propagator as well as the three- and four-gluon proper vertices using the ghost and gluon propagators 
obtained from the Gaussian wave functional [8] as input. A full self-consistent inclusion of the three- and four-gluon 
vertices will be the subject of future research. 

It is clear from the very beginning that eventually we have to truncate the tower of DSEs for the proper n-point 
vertex functions r„ as well as the equations of motion for the variational kernels 7„ following from the variational 
principle. For a systematic counting of the various diagrams we will assume a skeleton expansion. 

Previous variational calculations (using Gaussian wave functionals) were restricted to two (overlapping) loops in 
the energy (if), resulting in a one- loop gap equation for the gluon propagator. Restriction to two overlapping loops 
in the energy results in a bare (zero-loop) three-gluon kernel 73 and in a vanishing four-gluon kernel, 74 = 0. To get 
a 74 7^ 0, one has to include up to three loops in the energy. To keep the calculation sufficiently simple, we will keep 
only those three (overlapping) loop terms in the energy containing three- or four-gluon kernels. This will result in a 
bare (zero-loop) four-gluon and a one-loop three-gluon vertex. 

The organization of the paper is as follows: in Sec. II we present the DSEs of the Hamiltonian approach first for a 
general field theory and afterwards for Yang-Mills theory in Coulomb gauge. The full static (equal-time) propagators 
of the Hamiltonian approach are expressed in terms of proper vertex functions in Sec. III. In Sec. IV we specify 
our Yang-Mills vacuum wave functional and derive the corresponding DSEs for the gluon and ghost proper n-point 
functions. By means of these DSEs, the vacuum expectation value of the Hamiltonian is expressed in Sec. V in terms 
of the variational kernels of the vacuum wave functional. In Sec. VI these kernels are determined by minimizing the 
energy density. Finally, in Sec. VII we calculate the three- and four-gluon proper vertices using as input the ghost 
and gluon propagators from the variational calculations with a Gaussian wave functional. A short summary and our 
conclusions are given in Sec. VIII. 

II. DYSON-SCHWINGER EQUATIONS OF THE HAMILTONIAN APPROACH TO YANG-MILLS 

THEORY IN COULOMB GAUGE 



Consider a quantum field theory comprised of a collection of fields <f> — {4>i, 4>2, ■ ■ ■) and let \ip) be the exact vacuum 
state. All static (time-independent) Green's functions, i.e., vacuum expectation values (<j><j)...), can be calculated 
from the generating functional 



where j = (ji , fa , . . . ) stands for the collection of sources corresponding to the fields and we use the abbreviation 
3 " 4> = ji'Pi + .7202 + • • • • In the "coordinate" representation of the vacuum state (<j>\ip) = il>[<t>], the scalar product in 
(1) is defined by the functional integral over time-independent fields 4>(x) 



Furthermore, the integral in the exponent is over spatial coordinates x of the static fields </>(x). Expressing the vacuum 
wave functional in the form 1 



A. General DSE formalism of the Hamiltonian approach to quantum field theory 



Z[j] = <V>|e/^|V>> 



(1) 




(2) 




(3) 



1 As long as we ignore the 0-vacuum of Yang-Mills theory, the vacuum wave functional can be chosen to be real, which we will assume 
in the present paper. 
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the generating functional of the Hamiltonian approach to quantum field theory becomes 

Z[j] = Jv^e-SM+f'-*, (4) 

which is a standard generating functional of the d — 3-dimcnsional Euclidean quantum field theory defined by an 
"action" S[4>]. Here, this action is defined by the vacuum wave functional ip[4>] and will, in general, be non-local and 
non-linear. We therefore perform a Taylor expansion of the action functional S[<f>] in powers of the time-independent 
fields 0(x). The constant part 5[0] is fixed by the normalization of the wave functional and the linear part can be 
absorbed into the external source. It is then sufficient to consider expansions of S[<j>] starting at second order 



S[4] = \ J 12 4? + ^ J 73 4? + ■ 



(5) 



Restricting the expansion to second order yields a Gaussian wave functional (3) for which the functional integral in 
Eq. (2) can be explicitly carried out. This corresponds to the so-called mean-field approximation, where all higher 
order Green's functions of the field 0(x) are given in terms of the propagator {4><p). 

In many cases the mean-field approximation is, however, not sufficient. Going beyond the mean-field approximation, 
the functional integral in (2) can no longer be explicitly performed. However, we can calculate the desired Green 
functions by exploiting Dyson-Schwinger equation techniques. Starting from the identity 
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D0l( e -Wi*)=o (6) 



we can derive, in the standard fashion, a set of Dyson-Schwinger equations (DSEs) for the Green functions (4><f>- ■ •). 
This infinite tower of equations has to be truncated to get a closed system of equations, and further simplifying 
assumptions on the form of the interaction kernels 7„ entering the ansatz for the vacuum wave functional, see Eqs. (3) 
and (5), will be required. Nevertheless, this approach allows us to go beyond Gaussian wave functionals and calculate 
the static Green functions (4>4>- • •) in terms of the kernels 7„. By means of these static Green functions, the vacuum 
expectation value of the Hamiltonian (ip\H\ip) is expressed in terms of the kernels j n , which are then found by 
minimizing the energy density. 

In the Hamiltonian approach to quantum field theory one is not primarily interested in the generating functional 
(2) itself but in expectation values of observables, in particular of the Hamiltonian. For this purpose it turns out to 
be more convenient to generalize Eq. (6) to 

J V4> ^ (c-^J Kft]) = J Vcp A {r ft] K [4,]iP[if,]) = 0, (7) 
where Kft] is an arbitrary functional of the underlying field (p. 



B. Derivation of the DSEs for the Hamiltonian approach to Yang— Mills theory in Coulomb gauge 

Below, we apply the general Dyson-Schwinger approach to the Hamiltonian formulation of quantum field theory 
outlined above to Yang-Mills theory in Coulomb gauge (which also assumes Weyl gauge Aq = 0). Implementing the 
Coulomb gauge by the Faddeev-Popov method, the expectation value of a functional K[A] of the (spatial components 
of the) gauge field A is given by 

(K[A])= f VAJ[A]\m\ 2 K[A\- (8) 

Here, tp[A] = (A\if>) denotes the Yang-Mills vacuum wave functional restricted to transverse fields, diAf — 0, and 
J\A\ = Dct(G^ 4 1 ) is the Faddeev-Popov determinant with 

G- A x = (S ab dl - g Af (x)9f )<5(x - y) (9) 

being the Faddeev-Popov operator. Since we work only with spatial vectors, we will use only Lorentz subscripts. 
Furthermore, g is the coupling constant, A ab = f acb A c is the gauge field in the adjoint representation of the colour 
group, and f acb are the structure constants of the su(N c ) algebra. The functional integration in Eq. (8) runs over 
transverse field configurations and is restricted to the first Gribov region Q or, more precisely, to the fundamental 
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modular region [23]. Moreover, we assume that the wave functional ip[A] is properly normalized, (ip\ip) = (1) = 1. 
Writing the vacuum wave functional as in Eq. (3) 

|^[A]| 2 =:c- s ^l (10) 

and choosing 

K[A] =eJ> A , (11) 

Eq. (8) becomes the generating functional of the static Green functions of the (transverse) gauge field A. In the 
following, it will be convenient not to fix K [A] to the form (11) but rather to let K[A] be an arbitrary functional of 
the gauge field. Furthermore, to simplify the bookkeeping we will use the compact notation 



AH ( x i) = - A - B = A(l) 5(1) = J d d x A«(x) B?(x) , 



(12) 



such that a repeated label means summation over the discrete colour and Lorentz indices along with integration over 
the spatial coordinates. 

Consider now the following identity 

= L VA sXtT) { J[A] C ~ S[A] K[A] }> (13) 

which holds due to the fact that the Faddeev-Popov determinant J\A\ vanishes on the Gribov horizon <9fi, tacitly 
assuming that the considered functional K[A] does not spoil the vanishing of J^AJif [A] on <9f2. Eq. (13) with K[A) 
given by Eq. (11) becomes the ordinary DSE of the usual (Lagrangian based) functional integral formulation of Yang- 
Mills theory in Coulomb gauge [24] when the functional integration is extended over time-dependent gauge fields 
Apfci) and S[A] is chosen as the usual classical action of Yang-Mills theory. 
Working out the functional derivative in Eq. (13) yields the following identity 



5 In J SS[A] 



_SA(1) SA(1) 
The derivative of In J can be written as 

<51nJ S 



SA(1) 5A{1) 
where we have introduced the bare ghost-gluon vertex 2 



-TrlnG^ 1 =r (l;3,2)G^(2,3), (15) 



^(2,3) 



r o( l;2,3) = -^. (16) 



With this result, Eq. (14) can be cast in the form 



/ 5 JB K[A] \ = / 6 JM\ + f o( l ; 3, 2) (G A (2, 3) K[A]) , (17) 



which is the basis of the gluon DSEs, exploited below in the evaluation of (H). 
Introducing ghost fields in the usual way 

J\A\ = Dct(G A 1 ) = J VcVce~ £G ^ c , 

the expectation value (8) explicitly reads 

(K[A]) = J VcVcK[A] e s[A]-cG^c ; 



(18) 



(19) 



The bare ghost-gluon vertex To defined by Eq. (16) differs from the one of Ref. [4] by an overall sign. 
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and Eq. (17) can be written as 



' A] K [ A }) = ( 6 -^) +r (l;3,2) (c(2)c(3)K[A}) . (20) 



The bare vertex F is the lowest-order perturbative contribution [7] to the full ghost-gluon vertex T defined by 

(A(l) G A (2, 3)) = (A(l) c(2) 5(3)) = -D(l, l') G(2, 2') f(l'; 2', 3') G(3', 3) , (21) 

where 

D(1,2) = (A(1)A(2)) (22) 

is the gluon propagator and 

G(l,2) := (G A (1,2)> = (c(l)c(2)) . (23) 

is the ghost propagator. 

Eq. (20) (or equivalently Eq. (17)) is the basic DSE of the Hamiltonian formulation of Yang-Mills theory in Coulomb 
gauge, and we will refer to it as "Hamiltonian DSE" . Below we will exploit this equation to express the various static 
(equal-time) correlators occurring in the vacuum expectation value of the Hamilton operator by the variational kernels 
7„, Eq. (5), of the wave functional ip[A]. This requires appropriate choices of the so far arbitrary functional K[A\. 

III. EXPRESSING STATIC CORRELATORS THROUGH PROPAGATORS AND PROPER VERTEX 

FUNCTIONS 

Choosing the functional K [A] in Eq. (20) as 

K[A] =exp{j • A + c-rj + fj-c}, (24) 

where j and fj, r] are the gluon and ghost sources, we obtain the generating functional of the full static (equal-time) 
Green functions 

Z[j, V , fj] = (exp{jA + c V + fjc}) =: , (25) 

where W[j,fj,r]] is the generating functional of the connected Green functions 



5W 



s 3 w 



M , „ s 2 w 

^ = w »= ' mw> 



j=7)=T)=0 



(26) 



Sj(l)Sfj(2)5r](3) 
Introducing the classical fields as 3 



= - (A(l)c(2)c(3)) , etc. 



A _$W -__5W _SW 
Sj ' drj 5fj ' 

we can define the effective action T[A, c, c] through the Legendre transform 

r[A, c, c] + W\j, ri,fj]=j-A + c-ri + fj-c, (28) 

where the sources have to be expressed by Eqs. (27) in terms of the classical fields A, c, c. From the effective action 
Eq. (28), the sources are obtained as 

■ <*r ST ST 

i = sa> ' = Tc> * = -Tc- (29) 



3 With a slight abuse of notation, we employ the same symbol for both the classical fields and the quantum fields which are integrated 
over. No confusion should arise, since they never appear together. Furthermore, derivatives with respect to Grassmann fields are always 
left derivatives. 
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Using Eqs. (27), differentiation with respect to the gluonic source can be expressed as 

6 _ SA(2) S Sc(2) 5 Sc(2) 5 
5j(l) ~ Sj(l) 8A{2) + Sj(l) 5c(2) + Sj(l) Sc(2) 

8 2 W S 5 2 W 6 S 2 W 

+ 



(30) 



5j(l)5j(2) 5A{2) 5j(l)5fj(2) Sc(2) 8j{l)5r,{2) 5c{2) 



Similar expressions can be written for the derivatives with respect to the ghost sources. Differentiating Eqs. (29) and 
using Eqs. (27), (30) we can link, in the usual way, the connected Green functions (derivatives of W) with the proper 
vertex functions (derivatives of T). As an example, we explicitly show how to express the full ghost-gluon vertex 
T(l; 2, 3), defined in Eq. (21), by derivatives of the effective action. We start from the identity 



6(1,2) 



8T 



8 2 W 



S 2 T 



+ 



6 2 W 



S 2 T 



8 2 W 



8 2 r 



6r](l) 8c(2) 6r](l)6j{2') 8A(2')Sc(2) 6r](l)6f}(2') 6c(2')Sc(2) 6t](l)6r](2') Sc{2')Sc(2) ' 



(31) 



which can be derived along the line of Eq. (30). Differentiating Eq. (31) with respect to a gluonic source j(3) and 
using Eq. (30) yields 







S 3 W 



S 2 T 



6 2 W 



6 2 W 



S 3 T 



(5j'(3) < 577(l) ( 577(2 / ) 6c(2')Sc(2) 5r)(l)Sf)(2>) Sj{3)Sj(3') 5A(3')5c(2')5c(2) 
where the omitted terms vanish when the sources are set to zero. By means of Eq. (26) and of 

6 2 T[A,c,c] 



+ 



Sc(2)5c(l) 



G~\l,2), 



A=c=c=0 



the last relation can be expressed as 



= -G-\2, 2') (A(3)c(l)c(2')> + G(2', 1)D(3, 3') 



S 3 T[A,, 



6c(2') Sc(2) SA(3') 



A=c=c=0 



Comparison with Eq. (21) shows 



T(l;2,3) 



5 3 T[A, c, c] 



5c(3)8c(2)SA(l) A=e=c=Q 
Similarly, defining the n-gluon proper vertex function by 

S n T[A, c,c 



r n = r(i,2,...,n) = 



5A(1) 5A(2) ■ ■ ■ 5A(n) 



(32) 



(33) 



(34) 



(35) 



(36) 



A=c=c=0 



the full gluon n = 2, 3, 4, 5-point functions defined by Eq. (8) with K [A] = AA . . . can be expressed through the proper 
vertex functions as 



D(l,2)^(A(l)A(2))=T(l,2)-\ 

(A(l) A(2) A(3)) - -r(l', 2', 3') D(l', 1) D{2' , 2) D(3\ 3) , 

(A(l) A(2) A(3) A(4)) = D(l, 2)D(3, 4) + D(l, 3)£>(2, 4) + D(l, 4)D(2, 3) 
+ D(l', 1) D(2', 2) D(3', 3) D(4', 4)|-r(l', 2', 3', 4') 

+ D(5, 5') fr(l', 2', 5)r(5', 3', 4') + T(l', 3', 5)r(5', 2', 4') + T(l', 4', 5)r(5', 2', 3') 



(A(l) A(2) A(3) A(4) A(5)) = -Y(l' , . . . , 5') D(l, I')... D(5, 5') + 

+ [(A(1)A(2)A(6)) T(6, 3', 4', 5')D(3, 3')D(4, 4')D(5, 5') + 9 combinations] 

- [D(l, 1') (A(1')A(6)A(7)) T(6, 6')r(7, 7') (A(6')A(2)A(3)) (A(7')A(A)A(5)) 

+ 14 combinations] 

- [£>(1, 2) (A(3) A(4) A(5)) + 9 combinations] , 



(37) 
(38) 

(39) 



(40) 





FIG. 1: Expression of the gluon three-point function, Eq. (38), by means of three-point proper vertex function and propagators. 
Here and in the following, fat shaded gray blobs represent full Green's functions, small filled dots connected Green's functions, 
and small empty dots stand for proper vertex functions. 





FIG. 2: Expression of the full gluon four-point function, Eq. (39). The prefactors indicate the number of possible permutations. 



The proper n-point gluonic functions T(l,2, ...,n) (36) are by definition invariant with respect to a permutation of 
external legs, i.e., of the entries 1,2, ... ,n. Eqs. (38)-(40) are represented in diagrammatic form in Figs. 1-3. The 
prefactors in Figs. 2, 3 indicate the number of the possible combinations. Consider, e.g., the second diagram on the 
right-hand side of Fig. 3: out of the five external legs one can form ( 2 ) = 10 pairs of external legs attached to the 
right vertex. The remaining three external legs have to be attached then to the left vertex and thus do not add 
more possible combinations. In the third diagram there are five possibilities to select the external leg attached to the 
internal lines. From the remaining 4 external legs there are ( 2 ) = 6 possibilities to choose the two external legs at the 
right external vertex. This fixes also the external legs at the left vertex. The symmetry of the diagram with respect 
to the interchange of the two external vertices introduces an extra factor ^. Therefore this diagram occurs with the 
multiplicity 5 • 6 • 5 = 15. 

In a similar fashion one finds for the expectation value of two gauge fields and a ghost and anti-ghost field 



(A(l) A{2) c(3) c(4)> = (A(l) A(2) G A (S, 4)) = D(l, 2) G(3, 4) 
+ D(l', 1) D(2', 2) G(3, 3') G(4', 4) j-f (l', 2'; 3', 4') 
+ r(l',2',5)D(5,5 , )r(5 , ;3',4') 

+ f (l' ; 3', 5)G(5, 5')f (2'; 5', 4') + f (2'; 3', 5)G(5, 5')f (1'; 5', 4') 
where the two-gluon-two-ghost vertex is defined by 

6 4 T[A, c,c] 



F(l,2;3,4) = 



<Sc(4) <Sc(3) 5A(2)5A(1) 



A=c=c=0 



The diagrammatic representation of Eq. (41) is shown in Fig. 4. 

The last four-point function we need for the evaluation of (H) is the ghost four-point function 



(41) 



(42) 



(c(l) c(2) c(3) c(4)) = ( [G A (1, 2) G A (3, 4) - G A (1, 4) G A (3, 2)] ) , (43) 
which can be expressed in terms of propagators and proper vertices in the standard way, yielding 
(c(l) c(2) c(3) c(4)) = G(l, 2) G(3, 4) - G(l, 4) G(3, 2)+ 

+ G(l, 1') G(3, 3') f (5; l', 2')D(5, 5')f (5'; 3', 4') x 
x [G(2', 2) G(4', 4) - G(2', 4) G(4', 2)] + 

- f (!', 3', 2', 4') G(l, 1') G(2', 2) G(3, 3') G(4', 4) , 



(44) 




FIG. 3: Expression for the full gluon five-point Green function, Eq. (40). 




FIG. 4: Vacuum expectation value of two gauge fields and a ghost and an anti-ghost field after Eq. (41). 
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FIG. 5: Vacuum expectation value of two ghost and two anti-ghost fields, Eq. (44). 



where the four-ghost vertex is defined by 

f(l,3,2,4) = 



<s 4 r[A, c, c] 



<5c(4) 8c{2) <5c(3) $c(l) 
Eq. (44) is represented diagrammatically in Fig. 5 



(45) 

A=c=c=i) 



IV. THE VACUUM WAVE FUNCTIONAL AND CORRESPONDING DSES 

So far, all manipulations have been exact. In Sec. II B we have presented the Hamiltonian DSEs for arbitrary wave 
functionals. To proceed further, we have to make an ansatz for the form of the vacuum wave functional tp[A], which 
by Eq. (10) defines the "action" functional S[A]. 

In perturbation theory, the vacuum wave functional in the form (3), (5) has been determined up to order 0(g 2 ) 
by a solution of the Schrodinger equation, and the resulting expressions for the kernels 72, 73, and 74 are given in 
Ref. [25]. In the present non-perturbative approach, we will assume a wave functional of the form (3), (5) with an 
"action" functional to be given by 

S[A] = w(l, 2) A(l) A{2) + 1 7(1, 2, 3) A(l) A(2) A(3)+ 

3 - (46) 
+ i[7 (l,2,3,4)A(l)A(2)A(3)A(4). 

For historical reason, we have denoted ^7(1,2) by w(l,2). Consistent with our convention on the proper vertices 
Eq. (36), we will frequently use the shorthand -f n = 7(1, . . . ,n). The functions to = ^72, 73, and 74 are variational 
kernels which will be determined by minimization of the vacuum energy density. As discussed before, Eq. (46) can be 
considered as arising in leading orders of a systematic Taylor expansion of the "action" functional. 

Let us also stress that the ghost fields do not enter the Yang-Mills vacuum wave functional ip[A]. The ghost fields 
are auxiliary fields to represent the Faddecv Popov determinant in local action form. By the very definition of the 
ghost fields, Eq. (18), the ghost-gluon vertex in the "action" (i.e., the exponent of Eq. (19)) has to be the bare vertex 
and, in principle, there is absolutely no need to include ghost or ghost-gluon kernels as variational kernels in the 
wave functional. However, due to approximations to be introduced, for practical purposes, one might also include 
ghost or ghost-gluon vertices as variational kernels in the wave functional to improve the latter. If the exact gluon 
wave functional ip[A] were used, the variational principle would determine the ghost-gluon kernel as the bare one and 
higher ghost kernels to vanish. Therefore we will not include additional ghost vertices into the variational ansatz for 
the vacuum wave functional. 

By construction, the variational kernels j n (which are purely gluonic) are totally symmetric with respect to permu- 
tations of the overall indices. Furthermore, the wave functional defined by Eqs. (10), (46) is normalizable even when 
the restriction of the functional integration to the first Gribov region is ignored, provided the kernel 74 is positive 
definite, which we will assume for the moment and which later on will be confirmed by our calculations. With the 
action functional Eq. (46), the Hamiltonian DSE (17) becomes 

2w(l, 2) (A(2)K[A}) + \ 7(1, 2, 3) (A(2) A(3) K[A\) + 

+ 1 7(1, 2, 3, 4) (A(2) A(3) A(4) K[A\) = (^y ) + ?o(l; 3, 2) {G A {2, 3) K[A\) . 
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FIG. 6: Diagrammatic representation of Eq. (48). The empty square boxes denote the variational kernels 7„. 

*0* = W (1,2) ^ = 7(1,2,3) ^=7(1,2,3,4) 

FIG. 7: Variational kernels occurring in the exponent of the wave functional. 

Except for the non-locality of the variational kernels u>, 73, and 74, the functional Eq. (46) has the same structure 
as the ordinary Yang-Mills action. Therefore, the DSEs resulting from Eq. (47) will have the same structure as the 
DSEs of ordinary d = 3 Yang-Mills theory in Landau gauge, however with bare vertices replaced by the non-local 
variational kernels u, 73, and 74. Eq. (47) is our fundamental DSE for the Hamiltonian approach to Yang-Mills 
theory in Coulomb gauge. 



A. DSEs of gluonic vertex functions 



The first DSE is obtained by setting K [A] = 1 in Eq. (47). Using (A) = and the expression Eq. (38) for the 
three-point function yields the identity 

= 1 7(1, 2, 3)D(2, 3) - f (1; 3, 2)G(2, 3) 

2 1 (48) 
- - 7(1,2, 3, 4) T(2', 3', 4') D(2, 2') D(3, 3') D(4, 4') , 

which is diagrammatically illustrated in Fig. 6. Eq. (48) is not really a dynamical equation but rather a constraint, 
which can be used to simplify tadpole terms in the evaluation of higher-order DSEs. It is also easy to see that in 
lowest order perturbation theory each term in Eq. (48) vanishes separately. 

The DSE for the gluon propagator follows from (47) by putting K[A] = A, yielding 

Ml, 3) (A(3) A(2)) + \ 7(1, 3, 4) (A(3) A{A) A(2)) + 

_ (49) 
+ - 7(1, 3, 4, 5) (A(3) A{A) A(5) A{2)) = i(l, 2) + r (l; 4, 3) (G A (3, 4) A{2)) , 

where we have introduced the abbreviation 

2) = 5 a ^ t klk2 (xi) 5(xi - x 2 ) (50) 

and tk 1 h 2 ( x i) = 5k 1 k 2 ~ ®kl®kl I ^xi 1S ^ ne transverse projector. By means of Eqs. (38), (39), the three- and four-point 
functions in Eq. (49) can be expressed through the proper vertex functions T n . By multiplying Eq. (49) by the inverse 
gluon propagator Eq. (37) and defining 

D(l,2)- 1 =r(l,2)=:2f2(l,2), (51) 

Eq. (49) can be cast in the form 

0(1, 2) = 2) - £(1, 2) + X (l, 2) + 2) - 0a(l, 2) + &(1, 2), (52) 
where we have introduced the following loop terms 

£(1, 2) - ± 7(1, 3, 4) D(3, 3') D(4, 4') T(3', 4', 2), (53a) 

X (l, 2) - I f (1; 3, 4) G(3', 3) G(4, 4') f (2; 4', 3'), (53b) 

0i (1, 2) = \ 7(1, 3, 4, 5) D(3, 3') D(4, 4') D(5, 5') D(6, 6') T(4', 5', 6) r(3', 6', 2), (53c) 

02(1, 2) = 3^ 7(1, 3, 4, 5) D(3, 3') D(4, 4') D(5, 5') T(3', 4', 5', 2), (53d) 

0t(l,2) = |7(l,2,3,4)D(3,4). (53e) 
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FIG. 8: Diagrammatic representation of the DSE for the gluon propagator, Eq. (52). 







FIG. 9: Diagrammatic representation of the DSE (55) for the three-gluon proper vertex function. The factor 2 in front of the 
ghost loop accounts for the two diagrams differing in the direction of the ghost line. 

Eq. (52) is represented diagrammatically in Fig. 8, and is recognized as the usual DSE for the gluon propagator of 
Landau gauge Yang-Mills theory [26], except for the replacement of the bare Yang-Mills vertices (defined by the 
Yang-Mills Lagrangian) by the variational kernels j n (defined by the ansatz (46) for the vacuum functional), which 
are represented by open square boxes, see Fig. 7. Note that the gluon loop £(1,2) Eq. (53a) disappears when the 
three-gluon kernel 73 is absent from the exponential Eq. (46) of the wave functional Eq. (10). For a Gaussian wave 
functional (73 = 74 = 0) only the ghost loop x(l, 2) Eq. (53b) survives from the loop terms in the DSE (52), Fig. 8. 
Choosing K[A] = A(2)A(3) in Eq. (47) yields the DSE for the three-gluon vertex, 



Ml, 4) (A(4)A(2)A(3)) + - 7(1, 4, 5) (A{4) A(5) A(2) A(3)) + 



3! 



7(1,4, 5, 6) (A(4) A(5) A(6) A(2) A(3)) = T (l; 5, 4) (G A (4, 5) A{2) A{3)) . 



(54) 



By means of Eqs. (38)-(41), the first two terms on the left-hand side and the right-hand side of Eq. (54) can be 
expressed in terms of proper vertex functions r„. The explicit evaluation of the five-point function is quite lengthy, 
and we quote only the result. Restricting ourselves to terms involving up to one loop (which is sufficient to obtain 
the energy (H) up to three overlapping loops, see the introduction), and chopping off the external propagators we 
eventually find from Eq. (54) the DSE for the proper three-point vertex function F(l, 2, 3) 

T(l, 2, 3) = 7 (1, 2, 3) + 7 (1, 4, 5) D(4, 4') D(5, 5') D(6, 6') T(2, 4', 6) T(3, 5', 6') 
- f (1; 4, 5) G(4', 4) G(5, 5') G(6', 6) [f (2; 6, 4') f (3; 5', 6') + 2 3] 



(55) 



--7(l,4,5) J D(4,4')i?(5,5')r(4',5',2,3) 
+ f (l;4,5)G(4',4)G(5,5')f(2,3;5',4') 
- l - [ 7 (1, 2, 4, 5) D(4, 4') D(5, 5') T(4', 5', 3) + 2 3] , 

which is represented diagrammatically in Fig. 9. 

Analogously, one can derive the DSE for the four-gluon vertex. Restricting ourselves again up to three loops in the 
energy (H), one finds just the "tree-level" expression 



r(l, 2, 3, 4) =7(1,2,3,4) + ... 



(56) 



Any loop contribution to T 4 generates at least four-loop terms in the energy, which are beyond the scope of the present 
paper. 
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FIG. 10: Leading order DSE for the four-gluon vertex. 



->-...cx 



FIG. 11: Diagrammatic representation of the ghost DSE (58). 



B. The DSEs for the ghost propagator and the ghost-gluon vertex 

Inverting the denning equation of the Faddeev-Popov operator (9) one finds the following identity [4] 

G A {1, 2) = G (l, 2) - G A (l, 4) A(3) f (3; 4, 5) G (5, 2) , (57) 

where Go(l, 2) = [(— 9 2 ) _1 ](1, 2) is the bare ghost propagator and T is the bare ghost-gluon vertex denned in Eq. (16). 
Taking the expectation value of Eq. (57) and using Eq. (21) yields the usual DSE for the ghost propagator [4] 

G(l, 2)- 1 = G (l, 2)- 1 - f (3; 1,4) G(4, 4') D(3, 3') f (3'; 4', 2) , (58) 
which in momentum space reads 

O-'V-J + .g-f*, p - q ,p) *MS. G ,p - , (59) 

and which is represented diagrammatically in Fig. 11. 

To derive the DSE for the ghost-gluon vertex there are two possibilities. The first one is to multiply Eq. (57) by 
the gauge field and to take the expectation value of the resulting expression. This leads to 

(A(l) G A (2, 3)) = - (A(l) A(4) G A (2, 5)) f (4; 5, 6) G (6, 3) . (60) 

The remaining expectation value can be expressed in terms of proper vertices by means of Eqs. (21) and (41). After 
chopping off the external propagators, this results in the DSE for the ghost-gluon vertex shown in Fig. 12. Eq. (60) is 
exact, i.e., not truncated, but not very convenient for the evaluation of the energy density. A more convenient form of 
the DSE for the ghost-gluon vertex is obtained by putting K[A] = G A in our general Hamiltonian DSE (47), thereby 
using the chain rule for the derivative of the ghost Green's function and using Eqs. (41), (44) to express vacuum 
expectation values through propagators and proper vertex functions. The resulting equation reads at one-loop level 

f (1; 2, 3) = f (l; 2, 3) + 7 (1, 4, 5) D(4, 4') D(5, 5') G(6, 6') f (4'; 2, 6) f (5'; 6', 3) 
+ f (1; 4, 5) G(4', 4) G(5, 5') D(6, 6') f (6, 2, 4') f (6'; 5', 3) 

- \ 7 (1, 4, 5) £>(4,4')I>(5, 5') r(4', 5'; 2, 3) 

+ f(l;4,5)G(4',4)G(5,5')r(2,5',3,4'). 

and is shown in Fig. 13. The main difference between Eqs. (60) and (61) is that while Eq. (60) is exact, in Eq. (61) 
two-loops terms involving higher-order vertices are neglected. Nevertheless, to our purpose, calculating the energy 
up to three loops, Eq. (61) is more convenient. We will use this equation in Sec. V to simplify the expression for the 
kinetic and Coulomb energy. 




FIG. 12: Diagrammatic representation of the DSE (60) for the ghost-gluon vertex. 
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1 = 1+ + ~ 

FIG. 13: Diagrammatic representation of the DSE (61) for the ghost-gluon vertex. 

V. ENERGY DENSITY OF THE YANG-MILLS VACUUM 

The DSEs of the Hamiltonian approach derived in Sec. IV are not "equations of motion" in the usual sense, but 
rather connect the various Green functions with the kernels occurring in the ansatz for the wave functional, while 
these kernels themselves are at this point not yet fixed. Here is where the variational principle comes into play: we 
will now evaluate the expectation value of the Yang-Mills Hamiltonian for the wave functional (10) with the ansatz 
(46) and then minimize it with respect to the variational kernels j n . 

The Yang-Mills Hamiltonian in Coulomb gauge reads [27] 



h= J d d x [l l 7- 1 [A]n?(x).7[A]n?(x) + lj;-.(x)i^.(x) 

+ L. jd d x d d y J- 1 [A] p a (x) FJJ 6 (x, y) J [A] p b (y) . 



(62) 



Here, II£(x) = -i6/6A%(x) is the momentum operator, and Ffi = diA a - - d^Af + gj abc A\A c J , is the non-Abelian 
field strength tensor. The first two terms in Eq. (62) are the electric (kinetic) and magnetic parts of the ordinary 
Yang-Mills Hamiltonian restricted to the curvilinear "coordinate" space of Coulomb gauge. The third (Coulomb) 
term arises from the resolution of Gauss's law and describes the interaction of the non-Abelian colour charge (of the 
fluctuating gauge field) with density 

p a (x)=if(x)n, b (x) (63) 
through the non-Abelian Coulomb interaction kernel 

Ff (x,y) = [ { -bd)-\-d 2 ){-bd)-^ y . (64) 

The Yang-Mills Hamiltonian in Coulomb gauge Eq. (62) is a positive definite operator. Accordingly the energy 
(ip\H\ip) is bounded from below (by zero) and the variational principle is applicable. 
For later use, we rewrite the magnetic term of the Hamiltonian in the symmetrized form 

\ F?j = -\ A ^ A + | T ^ + ^ T ± A " > ( 65 ) 
where the interaction kernels are given in momentum space by 

T^ip, q, k) = i f abc [5 l3 (p - q) k + 5 jk (q - k) t + S kl (k - p)j] (66a) 

and 



T t fkl d ={f abe f Cde (Srk 




— Su Sjk) 


+r ce f bde (s ij 


Ski 


- Sjk Su) 


+f ade fb ce {6ij 


Sm 


- S ik Sji)} 



(66b) 

The bare four-gluon vertex, Eq. (66b), is independent of the momenta. 

A. Technicalities 

To evaluate the vacuum expectation values of the kinetic term and of the Coulomb Hamiltonian, it is convenient 
to perform an integration by parts in the gauge field. This leads to expressions of the form 
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where f[A] is a functional of the gauge field which does not contain any momentum operator. In principle, we could 
now explicitly write down the variations of the action (46) and evaluate the expectation value (67). Then one would 
recognize that some terms can be combined and simplified by using the DSEs (52), (55). Therefore, a more efficient 
way to evaluate the expectation value (67) is to use the Hamiltonian DSEs from the very beginning. 
Putting K[A] = 5S/5A(2)f[A] in the general DSE (17) we obtain 

*SW 6 S[ A ]f[A A = ,*m 



6A(1) 5A(2) Jl V \6A(1)5A(2) 

/5f[A] 6S[A]\ ~ /5S[A] " 

The last two terms on the r.h.s. of Eq. (68) can again be re-expressed through the DSEs (17), thereby putting 
K [A] = Sf[A]/SA and K [A] = G A , respectively, and using the definition of the bare ghost-gluon vertex r , Eq. (16). 
This results finally in the relation 



\<L4(1) <L4(2) /[ V \6A(l)5A(2) n 1 / \SA(1) SA(2) / 

+ fo(l;4,3)(^M G ,(3,4)) + f (2;4,3)(^ G ,(3,4)) ( 69 ) 

+ f (1; 4, 3) f (2; 6, 5) (f[A] [G A (3, 4) G A (5, 6) - G A (3, 6) G A (5, 4)] ) . 

We stress that this is an exact identity, which holds for f[A] being an arbitrary functional of the gauge field only, i.e., 
not containing the momentum operator. Furthermore, it will be sometimes convenient to express the last expectation 
value in terms of ghost fields 

(f[A] [G A (3, 4) G A (5, 6) - G A (3, 6) G A (5, 4)] ) = (f[A] c(3) c(4) c(5) c(6)> . (70) 

B. Kinetic energy 

After an integration by parts, the vacuum expectation value of the kinetic part of the Yang-Mills Hamiltonian (first 
term in Eq. (62)) can be expressed as 

F 1 f VA 7\A\ H[A] H[A] 1 l 5S[A] 5S[A] \ f71l 
E k --]VAJ[A] ^ - - ^ ^ . (71) 

The last expectation value has precisely the form of Eq. (69) with f[A] = 1 and the two external indices contracted. 
The terms in Eq. (69) involving functional derivatives of f[A] then vanish, and with the explicit form of the action 
(46) we find for the kinetic energy 



E k — 77 



2w(l,l) + 2^(l,l) + ro(l;4,3)ro(l;6,5)<c(3)c(4)c(5)c(6)) . (72) 



Here, <f>t is the gluon tadpole term occurring in the gluon DSE (52) and being defined by Eq. (53e). The ghost 
four-point function (cccc) occurring in the last term can be expressed by means of Eq. (44) in terms of propagators 
and proper functions. Contracting Eq. (44) with the two bare ghost-gluon vertices as in Eq. (72), one obtains 

f (1; 4, 3) f (1; 6, 5) (c(3) c(4) c(5) c(6)> = 4 X (1, 3)D(3, 4) x (4, 1) - [f (l; 5', 6') 

+ f (1; 4, 3)f (7; 3', 6')f (7'; 5', 4')G(3, 3')G(4', 4)D(7, 7') (73) 
+ f (l; 4, 3)P(3', 5', 4', 6')G(3, 3')G(4', 4)] G(6', 6) G(5, 5') f (l; 6, 5) . 

where x(l)2) is the ghost loop defined by Eq. (53b). The terms in the square brackets represent precisely the first, 
third, and fifth term of the right-hand side of the truncated DSE (61) for the ghost-gluon vertex. Therefore we can 
use Eq. (61) to rewrite the terms in the bracket in Eq. (73) in more compact form 

E k = \[lo{1, 1) + <M1, 1) - X (l, 1) + 2 X (1, 2)D(2, 3) X (3, 1) + ry c (l, 1) - 772(1, 1)] , (74) 
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FIG. 14: Diagrammatic representation of the contributions (75) to the kinetic energy density. 

where we have introduced the abbreviations 

2 Vc (l, 2) = f (1; 3, 4)G(3', 3)G(4, 4')G(5, 5')f (6'; 5', 3')f (7'; 4', 5)D(6, 6')D{7, 7') l3 (6, 7, 2) , 

~ ~ (75) 

2,72(1, 2) = r (l; 3, 4)G(3', 3)G(4, 4')r(5', 6'; 4', 3')D(5, 5')D(6, 6') 73 (5, 6, 2) , 

see Fig. 14. Eq. (74) can be slightly rewritten by using the gluon DSE (52) to eliminate the variational kernel u> (or, 
more precisely, the sum ui + <j> t ) in favour of the inverse gluon propagator Eq. (51). Assuming furthermore that the 
various loop terms Eqs. (53) are colour diagonal, e.g. 

£#(k) = (k)£(k), etc. (76) 

which is guaranteed by global colour invariance, we can express the kinetic energy in momentum space as 

Et __ W-w-i) v J dp { W-riril' + m _ Mp) + Mp) _ 2lfc(p) } , m 

where the terms are defined in Eqs. (53c) and (53d), and rj c is given in Eq. (75). Furthermore, in Eq. (77) we 
have introduced the abbreviation 

d d p 

ap = 



(2tt) c 



and V is the spatial volume which arises as in Ref. [4]. 

The ghost and gluon loop contributions defined in Eqs. (53) read in momentum space 

x(p) = a^-iKd-i) / iq F ^ c(p; q p ' ~ q) f T"(-p-> q. p q)G(q)G( P - q ) , (78) 

In the above equations, T and T 3 are, respectively, the full ghost-gluon and three-gluon vertices defined by Eqs. (35), 
(36). We have also introduced here the contraction of two colour and Lorentz tensor structures through transverse 
projectors as 

A o B := A^(p, q, k) t u (p) t jm (cO t kn (k) Bt c j-P, -q, -k). (80) 

Furthermore, the term r]2, Eq. (75), has been discarded, since it gives rise exclusively to higher-order contributions 
with more than three loops in the energy, which are beyond our truncation scheme. 

C. Magnetic energy 

In the notation of Eq. (65) the magnetic energy is given by 

Eb - \ (F%) = -\ (Ad- A) + | T 3 (A3) + £ T 4 (A±) . (81) 

The first term on the right-hand side of Eq. (81) can be expressed by means of the gluon propagator (22), while the 
second term can be expressed through the proper three-point function T3, Eq. (38) 

|T 3 (A 3 ) = -|T 3 or 3 (AA) 3 . (82) 
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FIG. 15: Contributions from the three- and four-gluon vertices to the magnetic energy, from left to right: second term in 
Eq. (83), Eq. (84), Eq. (85) and Eq. (86). The filled diamonds stand for the bare vertices T 3 , T 4 (Eqs. (66)) occurring in the 
magnetic part of the Hamilton operator, Eq. (65). 



In momentum space, these two energy contributions read 



n(p) 8 



v r A A (r 3 or3)( P ,q,k) 

— / dfc - FT7 - T7T7 - T7 ^- r <5(p + q + k). 



fi(p) fi(q) O(k) 



(83) 



Let us turn now to the last term in Eq. (81), (A 4 }. The four-point function is expressed by means of Eq. (39) in 
terms of gluon propagators and proper vertex functions. The disconnected terms in Eq. (39), i.e., the products of two 
gluon propagators, when contracted with the bare four-gluon vertex T4 (66b) results in 



16 



fi(p)fi(q) 



(84) 



which is the usual gluon tadpole term, which occurs already when a Gaussian wave functional is used. Notice that, 
since the q integral does not depend on an external momentum, we can replace q^qj — > ^ij i n the integrand, and 
using 



d(d - 3) + 3 



(d-lf 



we can rewrite Eq. (84) as 



,N c {Nl-l) (d-lf 



16 



d 



7 



V dpdq 



0(p) O(q) 



(84') 



Besides this, we get from the last term in Eq. (81) by using Eq. (39) also a contribution containing the proper 
four-point vertex function T4 



g 2 V 



dp dq dk dl- 



(T 4 o r 4 ) 



i6-4iy 1 * si{p) o(q) n(k) n(€) 

and a contribution containing two three-gluon vertices 



(5(p + q + k + £). 



g 2 V 



d [pqkt] - 



t$* r^ e (-P, -q, p + q) r^jk + t, -k, -I) 



32^(p) O(q) O(k) 0,(1) 0(p + q) 



(27r) d (5(p + q + k + £). 



(85) 



(86) 



The Lorentz indices in Eq. (86) are supposed to be contracted by transverse projectors, which we have not explicitly 
written down in order to prevent the equation from getting cluttered. The various interaction contributions to the 
magnetic energy given by Eqs. (83)-(86) are shown in Fig. 15. Notice that the diagrams with a four-gluon vertex 
contain already three loops. 



D. Coulomb energy 



After an integration by parts, as in the case of the kinetic energy, the vacuum expectation value of the Coulomb 
Hamiltonian, last term in Eq. (62), can be expressed as 



»- jVAJ[A] jd d xd d y 



A ac (x) 



Sip[A] 



iM?(x) 



Fa fry) 



<(y) 



6ip[A] 
iSAjj) 



(87) 
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In order to exploit our compact notation, we rewrite the colour charge density (63) as 

p(l) = R{l;2,3)A(2) T ^, (88) 

where 

i?(l; 2, 3) = -R(l; 3, 2) = f a ^ a " 5 l2l3 5( Xi - x 2 ) <5( Xl - x 3 ). (89) 
In this notation, the Coulomb energy (87) reads 

E c = ^ R(l; 3, 4) J2(2; 5, 6) (f a (1, 2) A(3) A(5)^ . (90) 

The remaining expectation value can, in principle, be expressed in terms of (so far unknown) higher order vertex 
functions. The required manipulations are, however, quite involved, and further simplifications are needed for practical 
reasons. Here we will again restrict ourselves to terms containing up to three overlapping loops in the energy. With 
this approximation we can factorize the Coulomb kernel Fa (see also Ref. [4] ) as 

E c ~ ^ R(l: 3, 4) R(2; 5, 6) (F A (1, 2)) ^(3) A(5) |M . (9 1) 

The Coulomb propagator (Fa) is discussed later in Sec. VIC. The remaining expectation value has precisely the form 
(67) with f[A] = AA, and from Eq. (69) we obtain 

(il ^ ^ = (saw A{3) A{5) ) + ( saSw) ^ A ™ 
+ Fo(4; 8 ' 7) (sW) [A{3) Am Ga{7, 8) ) + Fo(6; 8 ' 7) (^i) [j4(3) Am Ga{7 ' 8) ) (92) 

+ f (4; 8, 7) f (6; 8', 7') (A(3) A{5) [G A (7, 8) G A (7', 8') - G A (7, 8') G A (7', 8)] ) . 

When inserted into Eq. (91), this expression can be simplified by exploiting the colour antisymmetry of the vertices 
R, see Eq. (89): For this reason, the part of Eq. (92) symmetric with respect to the interchange of the indices (3 o 4) 
or (5 6) vanishes. Furthermore, the two terms on the r.h.s. of Eq. (92) with a single (bare) ghost-gluon vertex Tq 
yield identical contributions when inserted into Eq. (91). 

The first four terms on the right-hand side of Eq. (92) can be straightforwardly evaluated as in the preceding 
sections. With the explicit form of the "action" (46), one obtains for these terms 

i(3, 6) i(4, 5) + 2 f (4; 8, 7) i(3, 6) (A(5) G A (7, 8)) + 2 w(4, 6) (A(3)A(5)) + 
+ 7 (4, 6, 7) (A(7) A(3) A(5)) + \ 7 (4, 6, 7, 8) (A(7) A(8) A(3) A(5)) = 

- t(3, 6) [t(4, 5) - 4 X (4, 7)D(7, 5)] + 2w(4, 6)D(3, 5)+ 

- 7 (4, 6, 7) D(5, 5') D(3, 3') D(5, 5') T(3', 5', 7') + 1 7 (4, 6, 7, 8) (A(7)A(8)A(3)A(5)) , 

where we have used Eqs. (21), (22), and the definition of the ghost loop Eq. (53b). We are still left with the four-point 
function in Eq. (93), and with the six-point function in the last term in Eq. (92). To work out these terms, we notice 
that the Coulomb energy (91) can be diagrammatically represented as 




where the "blob" (including the four lines attached to it) represents the v.e.v. given by the last bracket in Eq. (91) and 
the double line stands for the Coulomb propagator (Fa)- Since we need the energy up to three loops, we should keep 
only those contributions to the "blob" which either factorize in two disconnected lines or which are either irreducible 
or at most one- particle reducible (no box diagrams). At this order, for the last term in Eq. (93) it is sufficient to keep 
from the expression for (A 4 ) given in Eq. (39) only the disconnected terms, yielding 

§7(4,6,7,8) (A(7)A(8)A(3)A(5)) = 

= \ 7 (4, 6, 7, 8) [D(3, 5) D(7, 8) + £>(3, 7) D(5, 8) + D(3, 8) D(5, 7) + ...} (94) 
= £>(3, 5) 2 «M4, 6) + 7 (4, 6, 7, 8) D(7, 3) £>(8, 5) + . . . 



17 



where we have used the definition of gluon tadpole 4>t (53e) and made use of the symmetry properties of the four-gluon 
kernel 74. The dots in Eq. (94) stand for terms which give rise to energy contributions with more than three loops. 
Following the same line of reasoning, the last term in Eq. (92) can be transformed to 

f (4; 8, 7) f 0(6; 8', 7') (A(3) A{5) [G A (7, 8) G A {7', 8') - G A {7, 8') G A (7', 8)] ) 
= f (4; 8, 7) f 0(6; 8', 7'){D(3, 5) ( [G A (7, 8) G A (7', 8') - G A (7, 8') G A (7', 8)] ) 

+ (A(3)G A (7',8'))(A(5)G A (7,8)) + ...} 
= D(3, 5) [-2 x (4, 6) + 4 X (4, 7)D(7, 8)x(8, 6) + . . .] 

+ 4D(3,7)x(7,6)D(5,8)x(8,4) + ... 

where we used Eqs. (21), (44) and discarded terms involving more than three loops in the energy. Collecting all 
terms given by Eqs. (92)-(95) and inserting the result into (91), we finally obtain the Coulomb energy to the desired 
(three-loop) order 

E c = ^i?(l;3,4)i?(2;5,6)F(l,2)x 

x {2 D(3, 5) [w(4, 6) + -M4, 6) - X (4, 6) + 2 x (4, 7) D(7, 8) X (8, 6) 

+ i(3, 6) [t(4, 5) - 4 X (4, 7)£>(7, 5)] + 4£>(3, 7) X (7, 6)D(5, 8) X (8, 4) ( 96 ) 
- 7 (4, 6, 7) D(5, 5') D(3, 3') D(5, 5') T(3', 5', 7') 

+ 7(4,6,7,8)D(7,3)D(8,5)}. 

where i^(l,2) = (F A {1,2)). For later use we rewrite this expression in momentum space. Exploiting the symmetries 
of the entries and expressing the variational kernel uj through the inverse gluon propagator fl, as we did for the kinetic 
energy, Eq. (96) can be cast into the form 



E c = 9 



2 N C (N* - 1) w t. _, [d - 2 + (p ■ q) 2 ] 



V JdpdqF(p + q) 



16 J 1 * ^ ^ fi(p)ft(q) 

x { [fi(p) - x(p) - fi(q) + X(q)] 2 + £(p) fi(p) + C(q) O(q)} 



/dp dg dl F(l) ^ (P + . q) /^/^ x (97) 

8 7 fi(p)0(q)0(p + q) 

x 7f™„ (P + q,^ - P, -q - I) T c 4(-P - q, P, q) 

+ v£ja P d q ae f^ ^ff / abc / ade 7^ - p, p, - q, i) 

where F(k) is the Fourier representation of the Coulomb propagator (-Fa), see Sec. VIC below. If we discard the 
three- and four-gluon kernels 73 and 74, which also removes the gluon loop £(p) Eq. (53a), this expression reduces to 
the Coulomb energy obtained in Refs. [4, 5] with a Gaussian wave functional, 



, iV c (iV c 2 - 1) T/ [ d - 2 + (p ■ q) 2 



B c = <T ' V /dpdg J F 1 (p + q 



16 J fi(p)n(q) (98) 

x{[0(p)-x(p)-0(q)+x(q)] 2 }, 

The new features introduced by the inclusion of the three- and four-gluon kernels will be studied in the subsequent 
sections. 



VI. DETERMINATION OF THE VARIATIONAL KERNELS 



In the previous section we have expressed the vacuum energy (ip\H \ ip) in terms of the variational kernels w, 73, and 
74 (occurring in our ansatz (10), (46) for the vacuum wave functional ip[A\) and of the proper vertices T n , T n . We 
now use the DSEs (52), (55), (56), and (61) to express the proper vertex functions r„, T n occurring in the energy in 
terms of the variational kernels 7„. We are then in a position to determine these kernels by minimizing (ip\H \ ip). To 
make the calculations feasible, we will resort to a skeleton expansion of (ip\H\ip), keeping at most three-loop terms. 
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As we will see, this is the minimum number of loops required to obtain a non-trivial four-gluon kernel 74. In the 
variation of the energy with respect to ^72 = u> and 73 we will restrict ourselves up to two loops terms in (ip\H\ip), 
which will be sufficient to get a non-trivial 73 and a one loop gap equation for u. 



A. Three- and four-gluon kernel 



Below we determine the three-gluon kernel in leading order in the number of loops. For this purpose, it is sufficient 
to keep up to two-loop terms in the energy. The relevant contributions come then from the gluon loop £ (79) occurring 
in the kinetic energy (77), and the magnetic energy contribution (83). These terms contain the three-gluon kernel 
73 either explicitly or implicitly via the three-point proper vertex T 3 , which by its DSE (55) is given in lowest order 
by the three-gluon kernel 73. All remaining terms of contain additional loops and will henceforth be discarded, 
resulting in the "tree- level" expression T 3 =73. Inserting this expression in Eqs. (79) and (83) and taking into account 
the symmetry of these kernels, the variation of these energy terms with respect to the three-gluon kernel 73 leads to 



_5_ 

<$73 



/ 



dpdqd£- 



<5(p + q + . 



'fi(p)fi(q)fi(£) 
which fixes the three-gluon kernel to 



(73 73) 



n( P ) + n(q) + n(£) 



- (73 °gT 3 ) 



= 0, 



%°/ fe c (p, q, k) n(p) + + 0(k) 



(99) 



(100) 



where the tensor T^(p, q, k) is defined in Eq. (66a). The obtained three-gluon kernel 73 is reminiscent of the 
perturbative one following from a solution of the Yang-Mills Schrodinger equation in leading order in the coupling 
constant g [25] 



(0)abc 



(p,q,k) = 



2glff fc e (p,q,k) 

IpI + |q| + l k l 



(101) 



except that the perturbative gluon energy |k| is replaced in Eq. (100) by the non-perturbative one Q(k). Note that, 
in principle, the Lorentz indices in Eqs. (100) and (101) are contracted with transverse projectors, which we did not 
explicitly write down, since they arise naturally as these kernels are always contracted with cither transverse gauge 
fields or the corresponding transverse Green functions. 

The variational determination of the four-gluon kernel 74 is technically somewhat more involved. The terms of (H) 
contributing to the variation with respect to 74 are those containing 74 either explicitly or implicitly via the DSEs 
(55), (56) for the proper three- and four-gluon vertices 1^, T4. Terms explicitly containing 74 are the <pi t 2 terms 
Eqs. (53c), (53d) of the kinetic energy Eq. (77), as well as the last term of the Coulomb energy (97). Terms containing 
T 3 or T4 are given by the gluon loop £ (79) in the kinetic energy (77) as well as by the magnetic energy contributions 
(83) and (85). All terms contributing to the variation of the energy with respect to 74 are collected below: 

^ 7(1, 2, 3, 4) D(l, 1') D(2, 2') £>(3, 3') 7 (l', 2', 3', 4) 

- ^ 7(1, 2, 3, 4) 7 (3', 4', 5) D(5, 5') 7 (5', 1, 2') D(2, 2') D(3, 3') D(4, 4')+ 



+ r(l,2,3)£>(l,l')£>(2,2') 



1 7 (1',2',3)-^(3,3')T(1',2',3') 



(102) 



- - 7 (1, 2, 3, 4) D(l, 1') D(2, 2') D(3, 3') D(4, 4') 7 (1', 2', 3', 4') 
+ ^ F(l, 2) R(l; 3, 4) R(2; 5, 6) 7 (4, 6, 7, 8) D(7, 3) D(8, 5) . 



For simplicity, we have not explicitly symmetrized this expression with respect to a permutation of the indices of 74. 4 
In the above expression we have used the DSE (56) to replace the four-gluon vertex T4 by the four-gluon kernel 74. 



4 Recall that, by definition, 74 is totally symmetric with respect to a permutation of indices. Strictly speaking, one should first sym- 
metrize Eq. (102) and then take the variation. However, the same result is more conveniently obtained by taking the variation of the 
unsymmctrized expression and symmetrizing afterwards. 
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To be consistent, we have to use the DSE for the three-gluon vertex, Eq. (55), to express the vertex function r 3 in 
(102) by the variational kernels, where it is sufficient to retain only terms involving 74 or r 4 , and the latter is to 
be replaced by 74 due to the DSE (56). Taking into account the symmetry properties of the quantities involved, by 
Eq. (55) we are led to make the following replacement in Eq. (102) 

T(l, 2, 3) -+ - \ 7(3, 4, 5) £>(4, 4') D(5, 5') 7 (4', 5', 1, 2) 

-7(3,1,4,5)D(4,4')D(5,5')7(4 , ,5',2). 

With this replacement, the variation of Eq. (102) with respect to 74 is now straightforward and yields after proper 
symmetrization with respect to external indices of 74 

[0(k 1 ) + 0(k 2 ) + n(k 3 ) + 0(k 4 )] 7$jftk 1 ,k 2 ,k 3 ,k 4 ) = 2ff 2 ^ b fe f 



7 ^(k 1 ,k 2 ,-k 1 - k 2 )t mn (k! +k 2 ) 7 ^(k 3) k 4) k 1 + k 2 ) 



(104) 



Hjri 

+ 7ffcm(ki, k 3 , -ki - k 3 ) t mn (ki + k 3 ) 7^„ e (k 2 , k 4 , ki + k 3 ) 

+ ^it (ki , k 4 , k, - k 4 ) t mn (ki + k 4 )7^ (k 2 , k 3 , ki + k 4 ) I 

_ 2g2 ^ r be f cde S . jSki [ 0(ki) _ 0( k2 )] F ( kl +k 2 )[0(k 3 ) - n(k 4 )] 

+ f ace f hde s lk s 3l [n(ki) - o(k 3 )] F(ki + k 3 ) [n(k 2 ) - n(k 4 )] 

+ f ade f bce SuS jk [fi(ki) - n(k4)] F(ki + k 4 ) [fi(k 2 ) - 0(k 3 )] J. 

F(p) is again the Fourier representation of the Coulomb propagator (Fa)- Eq. (104) yields a four gluon kernel 74, 
which is reminiscent of the perturbative one [25] except that the perturbative propagators and vertices are replaced 
by the full ones. 



B. Gap equation 

Given the explicit form of the energy functional (77), (83)-(86), (97), it is more convenient to use the DSE (52) 
for the gluon propagator to express the kernel w(p) in terms of the gluon energy f2(p), and vary the energy density 
with respect to fi _1 (p). The vacuum energy is given by closed loop diagrams, and the variation with respect to the 
gluon propagator Q _1 (p) reduces the number of loops by one. If the (gap) equation for O(p) is to be calculated up 
to one loop, it is sufficient to keep up to two loops in the energy. At this order, the energy density e defined by 
(H) =: V(d - l)(iV c 2 - l)e is given by 



-II 



dp P^_+Mp)- X (p)} 



9 2 N C f 73 
apaq 



n( P ) 

fi(p) + fi(q) + n(p + q 



! 73 - gT 3 



3!8(d-l)y ^ y Q(p)0(q)Q(p + q) 



(105) 



+ 16(d- C l) J SP ^ [d 2 + * ' 4)2]F(P + q) ' 



n(p)fi(q) 



see Eqs. (77), (83) and (98). In Eq. (105) we have discarded the tadpole term, Eq. (84) or (84'), since it represents 
an irrelevant constant, which disappears after renormalization. Except for the gluon loop (second term), the energy 
density Eq. (105) was already obtained in Ref. [4, 5], where a Gaussian wave functional multiplied by J~ X / 2 [A\ was 
used. The gluon loop is lost when a Gaussian wave functional is used, for which Green's functions with an odd number 
of fields vanish. 

In principle, in Eq. (105) the three-gluon kernel 7 3 (100) obtained from the variational principle Se/Sj3 — depends 
on f2(p). However, since the energy density e (105) with 7 3 given by Eq. (100) is already stationary with respect to 
variations of 73, we can ignore the implicit fi(p) dependence of e via 73. Variation of e (105) with respect to f2 -1 (k) 
then yields the gap equation 

n( p )2 = p 2 + x(p )2 _ /g(p) + /c(p ) 5 (106) 
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where x(p) 1S the ghost loop (78), 



'g(p) 



1 



d d q 



1 



4(d - l)(iV 2 -I) J (2ny n(q) fi(p + q) 
is the gluon loop contribution, and 

g 2 N c f d d q 



73 



9T3 - 73 



fi(p) + n(p + q) 



(107a) 



Hp) 



2(d-l)J (2tt) 



[»(q)-x(q) + x(p)] 2 -»(p) 2 
O(q) 



(107b) 



arises from the Coulomb term. Except for the gluon loop ic:(k), Eq. (106) is the gap equation found already in [4]. 

If we insert the expression derived above in leading order for the three-gluon kernel 73, Eq. (100), into the gluon-loop 
contribution Eq. (107a), this takes the form 



Hp) 



d- 



£ / dq 

1 J (2n) d 



2ft(p) + fi(q) + fi(p + q) E(p,q) 



[n(p) + n(q) + o( P + q)] 2 n (p + q) ' 



(108) 



where the function E arises from the contraction of the Lorentz structure of two three-gluon vertices upon imposing 
momentum conservation 



T 3 oT 3 | k= 

and is given explicitly by 

S(p,q)= [l-(p-q) 2 ] 



-p-q 



=:4iV c (iV 2 -l)£(p,q). 



(d-l)(p 2 +q 2 ) 



(d - 2)p 2 q 2 + (p ■ q)- 
(p + q) 2 



(109) 



(110) 



To exhibit the UV-behaviour of the various loop terms, let us consider them in leading order in perturbation theory 
[7, 25], where F(p) = 1/p 2 and x(p) = 0- I n this order we find with 3-momentum cut-off A 



Hp) = 
Hp) = 



9 2 N C 
(4tt) 2 
9 2 N C 
(4tt) 2 



4 a2 22 2] A 2 
3 A +15 P ! V 



A 



15 



p 2 ln- 



A 2 



(111a) 
(111b) 



One observes that the quadratic divergence in Eq. (106) cancels, and the sum of the logarithmic ones is consistent 
with the result of Lagrangian-based perturbation theory in Coulomb gauge [24, 28] 



fi(p) 



1 + ^Hp)-Hp)]) + --- = i + j0^ 



In- 



+ ... 



(112) 



where we have used dimensional regularization with d — 3 — 2s. 

To estimate the size of the contribution of the gluon loop (107a) to the gap equation (106), we use the gluon energy 
f2 obtained with a Gaussian wave functional in Rcf. [8] as input. For practical purposes, we fit the fi(p) obtained in 
Rcf. [8] to the formula 



fi(p) 




(113) 



which yields m\ ~ 0.36cr 2 , c 2 ~ 1.0tr c . The numerical results of Rcf. [8] and the fit to Eq. (113) are shown in Fig. 16. 
The gluon loop (107a) is UV divergent. In principle, the UV-divergent part is removed by the renormalization of the 
gap equation (106), which can be done analogously to Refs. [12, 29]. We therefore calculate 



ttnew(p) 2 = ^old(p) 2 ~ {HP)l^old] - HP)div) ■ 



(114) 



where Io(p)div is the (known) perturbative divergent part of the gluon loop, which has been subtracted, and Eq. (113) 
has been used for fl id(j>). The (inverse) gluon propagator Q new (p) Eq. (114) is shown in Fig. 17 together with the 
one obtained previously [8] with a gaussian wave functional fi /d(p)- The mismatch in the UV is due to the anomalous 
dimension developed by Q new (p) and absent in il id{p). As seen in Fig. 17, significant correction from the gluon loop 
arises in the mid- momentum regime, a behaviour which was observed also in Landau gauge [30, 31]. This a posteriori 
supports the use of the Gaussian wave functional for the description of the infrared regime. 
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FIG. 16: Numerical results from Ref. [8] for the gluon energy fi(p) and fit to Eq. (107a). 
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FIG. 17: Comparison of the gluon propagator obtained in Ref. [8] with a Gaussian wave functional (dashed line) with the 
corrected one from Eq. (114) (full line), and the lattice data from Ref. [9]. 



C. The Coulomb form factor 



With the explicit expression of the three- and four-gluon kernels, 73 Eq. (100) and 74 Eq. (104), on hand, we are 
left with three coupled equations: Eq. (59) for the ghost propagator, Eq. (61) for the ghost-gluon vertex, and the 
gap equation (106) for the gluon propagator. The final piece which is missing for closing this set of equations is 
the non-Abelian colour Coulomb potential F(l,2) defined by the vacuum expectation value of the Coulomb kernel 
Eq. (64) 

4) = (F A {1, 4)) = (G A (1, 2) G - 1 (2, 3) G A (3, 4)) . (115) 

The quantity g 2 F(l, 2) directly relates to the heavy quark potential [32] and is hence a renormalization group invariant 
quantity. Before we come to the evaluation of (Fa), let us remark that in practical application it is not necessary to 
explicitly solve the DSE (61) for the ghost-gluon vertex. Rather it is sufficient to replace the full ghost-gluon vertex F 
Eq. (35) by the bare one To Eq. (16). This approximation is motivated by the "non-renormalization" theorem for this 
vertex [33]. Although this theorem was originally proven [33] and confirmed on the lattice [34, 35] for QCD in Landau 
gauge, the arguments carry over to the present case of Coulomb gauge. A perturbative evaluation of the ghost-gluon 
vertex in Coulomb gauge shows indeed that its quantum corrections are finite and independent of the scale [7, 25] . 

The vacuum expectation value (F A ) is commonly expressed in terms of the Coulomb form factor /. This quantity 
measures the deviation of (Fa) from the factorized form (G A ) G^ 1 (Ga) and is defined in momentum space by 



F(p) =:G(p)/(p)p 2 G(p), 



(116) 
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where G — (Ga) is the ghost propagator (23). By taking the vacuum expectation value of the operator identity 

F A = ^-(gG A ), (117) 

the Coulomb form factor /(p) can be related to the ghost form factor d(p) and from the ghost DSE (59) the following 
(approximate) integral equation for the Coulomb form factor is obtained [4, 36] 

/(p) = 1 + J dq 1 -^ )2 (p - q) 2 G 2 (p - q) f(p q) . (118) 

In the derivation of this equation no assumption on the form of the vacuum wave functional enters, so this equation 
remains also valid for non-Gaussian wave functional. 

With the approximation T, Eq. (35), — > To, Eq. (16), Eqs. (59), (106), and (118) form a closed set of coupled 
equations for the ghost propagator G(p), the gluon energy £l(p), arL d the Coulomb form factor /(p), whose solutions 
provide the variational solution of the Yang-Mills Schrodinger equation. These equations have to be renormalized, 
which can be done in exactly the same way as in Rcfs. [4, 12, 29]. The numerical solution of this set of equations can 
be carried out as in the case of the Gaussian wave functional [4, 8]and is subject to future work. 

As already mentioned before, the use of a Gaussian wave functional misses the contribution of the bare three-gluon 
vertex and thus the gluon loop in the gap equation. Fortunately, as we have seen in the previous subsection, the gluon 
loop is IR subleading compared to the (included) ghost loop, and thus irrelevant for the IR behaviour of the Green 
functions. Therefore the results of Refs. [4, 8] remain fully valid in the IR. In particular, the gluon propagator does 
not change in the IR. The gluon loop does, however, matter in the UV and is responsible for the anomalous dimension 
of the gluon propagator, which enters the running coupling constant. 



VII. THE THREE- AND FOUR-GLUON VERTEX 



As we have seen, the gap equation differs from the one obtained in [4] only by the gluon loop contribution, which 
gives sizeable corrections in the mid-momentum regime. In the present section, we investigate the three- and four- 
gluon vertices T 3; 4 by using the ghost and gluon propagators determined with a Gaussian wave functional [8] as input. 
We will not resort to the tree-level result T n = j n but rather solve the corresponding DSEs to one-loop order. 



A. Solution of the truncated three-gluon vertex DSE 

The DSE for the three gluon-vertex at one-loop order is given by Eq. (55). Assuming ghost dominance (see 
for example Ref. [6]), we keep only the ghost- loop (the third term on the right-hand side of Eq. (55) and Fig. 9). 
Furthermore, we replace the full ghost-gluon vertex by the bare one, see Sec. VIC. After extracting the colour 
structure, the truncated DSE for the three-gluon vertex becomes 

r yfc (p, q, k) = lijk (p, q, k) + ig 3 N c J d£ G(£) G(£ - p) G(£ + q)£ t £ 3 (£ - p)k , (119) 

where momentum conservation k + p + q = 0is implied and we have omitted longitudinal terms, which, by definition, 
cannot enter the proper n-point vertex functions (36) of the transverse (Coulomb) gauge field. Possible tensor 
decompositions of the three-gluon proper vertex are given in Refs. [37, 38]. Here, for sake of illustration, we confine 
ourselves to the form factor corresponding to the tensor structure of the bare three-gluon vertex, defined by 

r (0) r 

where is the perturbative vertex given in Eq. (101). Restricting the kinematic configuration to the case where 
two external momenta have the same magnitude, i.e., 

q 2 = p 2 =p 2 , q-p = cp 2 , (121) 

the form factor Eq. (120) depends only on two variables, the magnitude p of the two momenta and the cosine c of the 
angle between them. Contracting Eq. (119) with the bare three-gluon vertex Eq. (101), thereby using the variational 
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FIG. 18: Ghost form factor from [8] and fit to Eq. (124). 



kernel 73, Eq. (100), and the tensor structure T 3 given in Eq. (66a), yields the following equation for the form factor 

fsA Eq. (120) 



f3A(p 2 ,c) 



2p+ ^/2p 2 (l + c) 
20(p 2 ) + r!(2p 2 (l + c)) 

2 N C 2+y/2(l + c) 



* 4 p(l- c)(9 + 8c + c 2 ) 
where we have introduced the abbreviation 



J dl G{£) G(£ - p) G{£ + q) B(£, p, q) , 



(122) 



£ 2 [p 2 (l-c) - (l + 2c)(£-q-^-p)] + 
-(l + c)[(£- q ) 2 + (£-p) 2 ] 
+ {£ ■ ci)(£ ■ p) [c(l + c) + 2 + • q - £ • p) 



(123) 



The apparent singularity at c = 1 in the coefficient in front of the integral in Eq. (122) is cancelled by the numerator, 
and the whole term is regular in the collinear limit. Also the singularities at £ = 0, £ = p, and £ = — q are integrablc. 

For the actual calculation of f^A, Eq. (120), we use the numerical result for Sl(p) and G(p) obtained in Ref. [8] 
(with a Gaussian wave functional) as input. For f2(p) we use the parametrization Eq. (113) while the ghost form 
factor d(p) = p 2 G(p) is parametrized as 



d(x) = a. 



1 



+ 



x 2 ln(x 2 + c 2 ) ' 



2 P 

x = — 

(To 



(124) 



with <7 C being the Coulomb string tension. The fit to the data shown in Fig. 18 yields a ~ 5 and c ~ 4. 

The form (124) of the ghost form factor embodies also the non-perturbative anomalous dimension, which in turn 
guarantees the convergence of the integral in Eq. (122). 

The numerical results for /3a(p 2 , c) are shown in Fig. 19 for some values of the cosine of the relative angle, c. The 
logarithmic plot (left panel) shows that the curves for different c have the same power law behaviour in the infrared 
region, namely p -3 , in agreement with the IR analysis of the ghost loop carried out in Ref. [38], according to which 
the IR exponent of the form factor j^A should be three times the one of the ghost dressing function Eq. (124). 5 The 
linear plot (right panel) shows that the form factor approaches unity in the high-momentum regime and changes sign 
in the mid momentum regime. 



5 In Ref. [6], due to the use of a Gaussian wave functional, the bare term escaped, and only the ghost loop was calculated for a different 
solution of the DSEs [4], which leads to less IR singular ghost and gluon propagators. 
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FIG. 19: Form factor /3a(p 2 ,c) defined in Eq. (120) for two momenta of equal magnitude p plotted for various values of the 
relative angle c. 
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FIG. 20: Form factor f^A of the three-gluon vertex for orthogonal momenta and comparison to lattice data for the d — 3 
Landau-gauge vertex [39]. The momentum scale is arbitrary and has been adjusted to make the sign change occur at the same 
point. The lattice data are shown by courtesy of A. Maas. 



In Ref. [39] the form factor /3A Eq. (120) of the three-gluon vertex was evaluated on the lattice for d = 3 Yang-Mills 
theory in Landau gauge. The result is shown in Fig. 20 and compares well to our result in low-momentum regime. 
In particular, in both studies, the sign change of the form factor occurs roughly at the same momentum where the 
gluon propagator has its maximum. The d = 3 Yang-Mills theory in Landau gauge can be interpreted as using the 
wave functional 

(125) 

in the Hamiltonian approach to d = 3 + 1 Yang-Mills theory in Coulomb gauge. This wave functional can be considered 
to represent the strong coupling limit of the true vacuum wave functional [40]. In Ref. [41] it was shown by means 
of lattice calculations that this wave functional yields static propagators which in the IR compare well with those 
obtained in d = 3 + 1 Yang-Mills theory in Coulomb gauge. Thus we can conclude from Fig. 20 that our results 
compare favourably with the lattice data. 



ip[A] = Af exp 



B. Estimate of the four-gluon vertex 



Due to its DSE (56), the four-gluon vertex is given in leading order (neglecting loops) by the variational kernel 74 
Eq. (104) determined in Sec. VIA. A form factor for the four-gluon vertex can be introduced along the same line of 
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FIG. 21: Form factor /4^(p 2 ) of the four-gluon vertex at the symmetric point. The leading contribution p/Q(p) is also shown. 
Eq. (120) 

r (o) o 

hA := r (° 4 ) IrW ' (126) 

1 4 O 1 4 



where rf^ is the perturbative four-gluon vertex, which is obtained from Eq. (104) by the replacements fi(p) 
F(p) —¥ 1/p 2 . We consider here the form factor f 4 A at the symmetric point, where 



P? 



2 2 
Pi = P ■ 



P 4 ' P, 



-3P 2 (i^j). 



(127) 



For this kinematic configuration the terms in Eq. (104) stemming from the Coulomb interaction do not contribute to 
this vertex, and one finds for the form factor (126) 



Ua(p 2 ) 



p 171 - m0[g o + g(p 2 )} + 8704 g g(p 2 ) 



where 



n(p2) 



g(p 2 ) - 



171 - 1920 5 o + 8704 g 2 



2n( P 2 ) + na P 2 ) 



.90 = g{p 2 ) 



(128) 

(129a) 
(129b) 



The function / 4j 4 Eq. (128) is shown in Fig. 21. The function multiplying p/Q,(p 2 ) in Eq. (128) is of 0(1), and the 
form factor at the symmetric point Eq. (127) can be fairly well represented by the gluonic dressing function p/VL{p 2 ) 
alone. 

In the above calculation of the four-gluon vertex we have neglected all loop diagrams, in particular the ghost loop. 
From analogous investigations in Landau gauge [42] one may expect that the ghost loop dominates the IR behaviour 
of the four-gluon vertex as it did for the three-gluon vertex. We will defer this issue to future research. 



VIII. SUMMARY AND CONCLUSIONS 



We have presented a general method to treat non-Gaussian wave functionals in the Hamiltonian formulation of 
quantum field theory. By means of well-established Dyson-Schwinger equation techniques, the equal-time Green 
functions and, in particular, the expectation value of the Hamiltonian are expressed in terms of kernels occurring 
in the exponent of the vacuum wave functional. These kernels are then determined by the variational principle 
minimizing the vacuum energy density. The method was applied to Yang-Mills theory in Coulomb gauge using a 
vacuum wave functional which contains up to quartic terms in the exponent. In leading order (in the number of loops) 
the cubic and quartic interaction kernels obtained from the variational principle are reminiscent of the corresponding 
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expressions obtained in leading order of a perturbative solution of the Yang-Mills Schrodinger equation, except that 
the unperturbed gluon propagators are replaced by the full ones. We have estimated these interaction kernels 73.4 
and the corresponding proper gluon vertex functions ^.4 by using the gluon and ghost propagators found in the 
variational approach with a Gaussian wave functional. The resulting three-gluon vertex compares fairly well to the 
available lattice data obtained in d — 3 Landau gauge. The (gap) equation of motion obtained from the variation 
of the energy density with respect to the gluon propagator contains the gluon loop, which was missed in previous 
variational approaches due to the use of a Gaussian wave functional. We have shown that the gluon loop gives a 
substantial contribution in the mid-momentum regime while leaving the IR sector unchanged. Furthermore (together 
with the contribution from the non-Abelian Coulomb interaction already fully included in previous studies [4]), it also 
provides the correct asymptotic UV behaviour of the gluon propagator in accord with perturbation theory [7]. The 
approach developed in the present paper allows a systematic treatment of correlations in the Hamiltonian approach 
to interacting quantum field theories (analogous to the so-called "exponential 5"' method in many-body physics) 
and opens up a wide range of applications. In particular, it allows us to extend the variational approach from pure 
Yang-Mills theory to full QCD. This will be subject to future research. 
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